Revealing the role of liquid preordering in crystallisation of supercooled liquids

The recent discovery of non-classical crystal nucleation pathways has revealed the role of fluctuations in the liquid structural order, not considered in classical nucleation theory. On the other hand, classical crystal growth theory states that crystal growth is independent of interfacial energy, but this is questionable. Here we elucidate the role of liquid structural ordering in crystal nucleation and growth using computer simulations of supercooled liquids. We find that suppressing the crystal-like structural order in the supercooled liquid through a new order-killing strategy can reduce the crystallisation rate by several orders of magnitude. This indicates that crystal-like liquid preordering and the associated interfacial energy reduction play an essential role in nucleation and growth processes, forcing critical modifications of the classical crystal growth theory. Furthermore, we evaluate the importance of this additional factor for different types of liquids. These findings shed new light on the fundamental understanding of crystal growth kinetics.


Yuan-Chao Hu 1 & Hajime Tanaka 1,2
The recent discovery of non-classical crystal nucleation pathways has revealed the role of fluctuations in the liquid structural order, not considered in classical nucleation theory. On the other hand, classical crystal growth theory states that crystal growth is independent of interfacial energy, but this is questionable. Here we elucidate the role of liquid structural ordering in crystal nucleation and growth using computer simulations of supercooled liquids. We find that suppressing the crystal-like structural order in the supercooled liquid through a new order-killing strategy can reduce the crystallisation rate by several orders of magnitude. This indicates that crystal-like liquid preordering and the associated interfacial energy reduction play an essential role in nucleation and growth processes, forcing critical modifications of the classical crystal growth theory. Furthermore, we evaluate the importance of this additional factor for different types of liquids. These findings shed new light on the fundamental understanding of crystal growth kinetics.
Crystallisation is one of the most ubiquitous yet mysterious phase transitions in nature. Understanding its physical mechanism is of practical interest in many fields, such as atmospheric science, nanoscale electronics, protein engineering, drug production, and physical metallurgy. Crystallisation usually proceeds by nucleation and subsequent growth, as described most famously by the classical nucleation theory (CNT) 1 . According to CNT, the equilibrium crystalline phase is nucleated randomly in a homogeneous disordered liquid by density fluctuations. Once the size of nuclei exceeds the critical size, they grow. The key thermodynamic factors for crystal nucleation are the free energy gain upon crystallisation, i.e., the chemical potential difference between the liquid and crystal phase Δμ, and the free energy cost proportional to the interfacial tension γ associated with the formation of a new interface between the two phases. The balance of these two factors determines the critical nucleation free-energy barrier ΔG and the critical nucleus size R c .
Although CNT serves as an essential framework to understand crystallisation, the discovery of non-classical pathways for crystallisation challenges its general validity [2][3][4][5][6][7][8][9][10][11] . This is mainly because of the assumptions involved in the simple CNT. For example, it was argued that the crystal phase nucleated from the liquid is not necessarily the thermodynamically most stable one but can be the one having minimal free energy difference with the liquid phase (known as the Ostwald step rule) 12 or the one with the lowest free energy barrier of formation 13 . It was also claimed based on the Landau theory that the body-centred cubic (bcc) phase should be nucleated preferentially and then transformed to the stable phase in simple atomic liquids 14 . Such two-(or multiple-)step crystallisation has been widely reported from numerical simulations and experiments, but still with controversy 2,3,6,15 . We may regard these scenarios as approaches from the crystal side.
On the other hand, recent studies of structural properties of glassforming liquids 16 have revealed that the supercooled liquid state is no longer homogeneous as assumed by CNT, and crystal-like angular order spontaneously develops in a supercooled state for systems suffering from only weak frustration against crystallisation. Such structural fluctuations have been found to play a dominant role in crystal nucleation, for example, by unsupervised learning methods 17 . That is, crystallisation begins with the enhancement of crystal-like bond orientational order, and then translational (density) ordering follows 4,7,10,16,18,19 . Thus, crystal nuclei are not born randomly but induced in regions of high crystal-like bond-orientational order in a supercooled liquid. In other words, the crystalline phase prefers to nucleate from preordered regions with local orientational symmetry consistent with the crystal. We may say that there is essentially no homogeneous nucleation in a strict sense. Nucleation always occurs through the "continuous" ordering, starting from orientational and then positional one 16,19 .
So far, this scenario has been mainly confirmed for relatively simple liquids with pair-additive isotropic potentials, such as hard spheres. How generally it is valid for more complex liquids such as metallic alloys is one issue. However, a much more crucial question is the role of such liquid structural preordering in crystal growth. Recently, fast crystal growth far beyond the prediction of CNT has attracted considerable attention [20][21][22][23][24] . According to classical crystallisation theory, the crystal growth rate is determined only by the diffusion constant D and the driving force of crystallisation Δμ, independently of γ. Considering preordering of supercooled liquids, it is not clear at all whether this is true.
In this work, to address these fundamental issues, we carry out extensive molecular dynamics (MD) simulations with NiAl (Ni 50 Al 50 ) as a model metallic alloy. The benefits of NiAl are two-fold: (1) its equilibrium crystal structure is B2, which is a bcc (body-centred cubic)-like structure consisting of two simple cubic interpenetrating sublattices; (2) spontaneous crystallisation is observable in the computational timescale. We find that the supercooled metallic liquid exhibits a nonclassical crystallisation pathway, generalising the findings from hardsphere-like systems with pair-additive potentials. More importantly, we show that both crystal nucleation and growth are mediated by local bond orientational order fluctuation rather than density fluctuation. Notably, by developing a novel order-killing strategy (OKS), we can control the degree of structural ordering in a supercooled liquid and consequently manipulate the crystallisation kinetics. Surprisingly, this method allows us to tune not only crystal nucleation but also growth rates to a large extent. We further show that the influence of structural ordering on crystallisation is through tailoring the liquid-crystal interfacial energy. This finding unambiguously demonstrates the crucial role of liquid preordering in crystallising metallic systems (including both nucleation and growth). The finding that liquid preordering affects crystal growth indicates the importance of the interface tension in crystal growth, unlike the prediction of the classical theory. Thus, an essential modification to the classical theory of crystal growth is necessary. We phenomenologically add the interfacial energy-related factor for estimating the impact of structural ordering on the growth rate. Then, we evaluate this factor for eight different systems with different bondings and crystal structures and discuss the implication of our results from the structural order type developed in these liquids.

Structural orderings during crystallisation
We start by studying the spontaneous crystallisation process in NiAl by annealing the supercooled liquid at the nose temperature T nose of its Time-Temperature-Transformation (TTT) curve (see Methods, section 'MD simulations of NiAl crystallisation'). Since the nucleation event itself is stochastic ( Supplementary Fig. 1a), we randomly follow one trajectory but find similar results from independent simulations ( Supplementary Fig. 2). The local structural features are characterised by bond orientational order parameters (see Methods, section 'Structural characterisation'). Figure 1a shows the time dependence of the fraction of crystallised atoms, together with that of preordered atoms with bcc-like Bcc-like local orientational ordering is the preordering linked to the crystal in this system because the symmetry is compatible with the equilibrium crystal. Fcclike and hcp-like preorders are almost negligible. There are always more atoms with ico-like order than those with crystal-like orders before crystallisation happens, i.e., in the incubation period. The temporal evolution of number density ρ is also included, which almost perfectly coincides with that of the crystallised atoms. symmetry. The time evolutions of other structural orders, with fcc-like, hcp-like, and ico-like symmetries, are also shown for comparison (fcc: face-centred cubic, hcp: hexagonal close-packed, ico: icosahedral).
Here, bcc-like, fcc-like, and hcp-like orders are categorised as crystallike since their symmetries are crystalline. We can find that the bcc-like one is dominant among these crystal-like orderings, and the others are almost negligible. It was proposed based on the Landau-type theory 14 that the metastable precritical nuclei with bcc order are formed first, and then their cores transform to the stable crystalline (e.g., fcc) phase in the growth process. Such behaviour was indeed observed for a Lennard-Jones liquid 15 . This kinetic pathway is known as a non-classical (two-step) crystallisation pathway. In our case, bcc-like preorders formed in the supercooled liquid share the same symmetry with the stable B2 crystal phase, whereas the other crystal-like preorders (fcclike and hcp-like) are negligible. Whether the bcc preordering is due to the Alexander-McTague scenario 14 or preferred by the interaction potential is an interesting question. We think that the latter may be the case. Anyway, our observation excludes the formation of a metastable crystalline phase during crystallisation. In Fig. 1, we can also see a plenty of ico-like clusters. This ordering may be favoured by the entropic effect since the icosahedral configuration provides high local packing capability, i.e., high vibrational entropy 25 . However, there is no quasicrystal formation mediated by ico-like preorders in our system, suggesting the strong internal geometrical frustration. Furthermore, we can see the strong competition between crystal-like and ico-like orderings from their mutually exclusive spatial distributions. Thus, ico-like orders only tend to prevent crystallisation without forming quasicrystals.
We can see crystallisation starting from approximately 5.0 ns. There are no crystallised atoms in the incubation period (t ≤ 5.0 ns), but always about 5% of the total atoms fluctuate as bcc-like preorders in the supercooled liquid. This conveys crucial messages. Firstly, the supercooled metallic liquid is not homogeneous as assumed by CNT but with local bond orientational order fluctuations heterogeneously in space. Secondly, these preorders are transient and have a finite lifetime. Thus, they are intrinsic to the supercooled liquid state. These features are similar to hard-spherelike systems 4,18 but still differ in that we find only bcc-like preordering in the metallic system 7 . This result can be understood as a consequence of the fact that preorders tend to share the symmetry of equilibrium crystals.
The atomic-scale features during crystallisation, as illustrated in Fig. 1b-d, provide further intriguing information. There are tiny preordered clusters fluctuating in space in the supercooled liquid state before crystal formation. The crystal nucleus is preferentially born within these preordered regions, initiated by local bond orientational order fluctuation but neither randomly as assumed by CNT nor by density fluctuation as assumed by the density functional theory. This is also revealed in Fig. 1a, where we can directly compare the number density change and bcc-like preorder. We can see that there is already a prominent amount of bcc-like preorder in the liquid state before crystal nucleation, which starts to increase before the onset of the density increase. We can see similar results in two other independent simulations (see Supplementary Fig. 2). We also stress that the local density field cannot detect these preordered regions (see Supplementary Fig. 3). So, we may say that the homogeneous nucleation of crystals occurs "heterogeneously in space". More precisely, the crystal nucleation is mediated by crystal-like preorders spontaneously formed in a supercooled liquid state. Since the nucleated crystal phase is bcclike, the same as the thermodynamically most stable one, it might look consistent with the classical scenario. However, the observation described above shows that crystallisation still proceeds in a two-step manner, more precisely, through sequential angular-then-positional ordering 4,16,18,19 . This fact indicates the generality of this physical picture to crystal nucleation in various materials, although each system has a specific preorder favoured by the interparticle interaction and entropy.
More importantly, Fig. 1 shows that the crystal nucleus grows while surrounded by preordered structures due to the interface wetting effect. This indicates a possibility that liquid preordering mediates not only crystal nucleation but also crystal growth. Therefore, the properties of the liquid-crystal interface can also be crucial in determining the crystal growth rate. However, such a possibility has been overlooked so far, e.g., in the Wilson-Frenkel (WF) theory 1,26,27 . Thus, we also focus later on how the preordering affects the crystal growth kinetics.

Tuning crystal nucleation time
For unveiling the role of such liquid preordering in governing the crystallisation kinetics and elucidating the underlying physical mechanism, we develop the OKS method (see Methods, section 'Order killing strategy') to intentionally control the degree of preordering in the supercooled liquid state. In brief, we check the local bond orientational orders of the liquid at a time window of Δt MD and perform MD simulations iteratively with and without a harmonic bias potential. The bias aims to revert the targeted preorders to a disordered state, i.e., killing specific preorders. The shorter Δt MD more strongly suppresses the targeted preorders in a supercooled liquid. In this way, we can suppress local bond orientational order fluctuations in a controllable manner.
We first show how the nucleation time at T nose is controlled by the structural orderings tuned by Δt MD in Fig. 2. Notably, crystal-like preordering plays a decisive role in initiating crystal nucleation. In standard isothermal annealing of the chosen crystallisation trajectory, nucleation begins from about 2.0 ns. However, when the fluctuation of the bcc-like ordering is suppressed, nucleation can be delayed considerably. The liquid with a shorter Δt MD , that is, fewer bcc-like preorders, requires a longer time to nucleate the crystal phase. Remarkably, the time required for nucleation increases to 1078 ns when Δt MD is decreased to 0.2 ns. This incubation time is more than 500 times slower than the one in the normal condition. Further decrease in Δt MD would lead to even a longer incubation time. Now, we consider the roles of ico ordering, which is another important type of ordering in our NiAl system. Topological orderings like local ico order in three dimensions are thought to be crucial in preventing crystallisation and glass formation from metallic The time dependence of the fraction of the crystallised atoms is shown at different conditions. Crystallisation originally starts from 2.0 ns during isothermal annealing (Normal MD). By killing the preordering periodically every Δt MD using the bias potential, the nucleation time has been remarkably changed. The shorter Δt MD is, the longer time is needed for nucleation to start. In contrast, the suppression of local icosahedral ordering (Δt ICO MD = 1:0 ns, dashed line) shows little influence on the incubation time of crystal nucleation. All the biased simulations start from the same condition from Normal MD. liquids [28][29][30][31][32][33] , which can be dated back to Frank's seminal work 34 . In addition to the number density, the spatial extendability of ico order has been suggested to be critical 35,36 . The formation of medium-range ico order, i.e., connected icosahedra, tends to suppress the development of crystal-like preorders more significantly than isolated icosahedra. The presence of spatially more extended ico ordering indicates its stronger ico ordering tendency, which overcomes internal structural frustration. Noting that crystallisation is controlled by thermodynamic competition between crystal-like and ico-like orderings, this means that crystallisation suffers from stronger frustration. Thus, the crystallisation kinetics should depend on the degree of such ico-like preorder suppression. It was shown by Desgranges and Delhommelle 36 that this spatial icosahedral extendability may induce a nonmonotonic temperature dependence of the 'crystal' nucleation barrier in supercooled metallic liquids. We also mention another possible influence of ico ordering on the fate of a liquid. When there is very strong tendency of ico ordering, their spatially extended connection may lead to the formation of metastable quasi-crystals 28,37 . Under such a circumstance, the 'crystal' nucleation is disfavored (the crystal nucleation barrier increases); instead 'quasicrystal' nucleation is favoured (the quasicrystal nucleation barrier decreases). Therefore, the spatial extent of ico-like order generally has significant effects on crystallisation. A critical factor is whether the symmetry is compatible with the (meta)stable phase to be formed or not.
In NiAl liquids studied here, when we eliminate the ico-like ordering instead of bcc-like ordering (see Methods, section 'Order killing strategy'), the nucleation time is almost unchanged and comparable to the normal one (Fig. 2). Thus, we conclude that the suppression of crystal-like order by ico-like order is not significant; in other words, the crystal-like preorder primarily controls crystal nucleation (see Fig. 1c). This minor influence of periodically killing only ico-like orders on crystallisation is because crystals are nucleated exclusively in regions of high crystal-like preorders. This explains the negligible influence of the ico-like order killing on the crystal nucleation time. However, if there is much stronger tendency of ico-like ordering, it may affect the crystallisation. How significantly ico-like orderings affect the crystallisation kinetics should depend on how strongly they suppress the crystal-like orderings. This is an interesting topic for future study.
This result does not rule out the importance of icosahedra but rather helps unveil the mechanism of how icosahedra affect the crystallisation kinetics microscopically, as discussed above. Due to the incompatibility of local symmetry between icosahedra and crystal-like preorders, the appearance of icosahedra automatically suppresses and prevents crystal-like preordering. These two types of orders are mutually exclusive in a supercooled liquid. Therefore, the way by which icosahedra influences crystallisation is through suppressing preordering and changing the liquid-crystal interfacial energy 28,38 . We also note that in our simulations, both short-range and medium-range orders are considered. Because we focus on detecting the centres of crystal-like or ico-like preorders, we destroy these orders without differentiating whether they are isolated or connected. Separating the effects of short-range and medium-range orders is another interesting topic for future study.

Tuning crystal growth rate
Next, we investigate the role of preordering in crystal growth at T nose by OKS (see Methods, section 'Order killing strategy'). To study the rate of crystal growth, we first identify the critical nucleus size N 0 c in the normal state as 51 ± 15 by using the seeding method 38 . In the course of spontaneous crystallisation, configurations with various sizes of crystal seeds can be obtained. For example, we choose a seed of 375 atoms embedded in the supercooled liquid as the starting state. As shown in Fig. 3a, the number of the crystallised atoms N xtal increases abruptly during normal isothermal annealing. However, when the preordering is suppressed via OKS, the crystal growth slows down remarkably. A shorter Δt MD , i.e., more frequent order killing, leads to slower crystal growth. The crystal phase can even be destroyed if the preordering is strongly suppressed (see, e.g., the case of Δt MD shorter than 25 ps in Fig. 3a). Decreasing Δt MD reduces the crystal growth speed considerably without melting, thanks to the large size of the crystal seed. However, the seed can even be dissolved when the crystallike preordering is strongly suppressed by very short Δt MD . b, The time dependence of N xtal in the following two-step killing procedures. We take crystallisation states from the case of Δt MD = 32 ps shown in a. Then, we reduce Δt MD to 10 ps at different growth stages, resulting in the more frequent removal of crystal-like preorders in this process. Since the larger crystal nucleus accompanies the larger amount of preorders, a shorter Δt MD is required to suppress the growth efficiently. The crystal growth rate can be effectively tuned by the degree of preordering through Δt MD . c, The dependence of the normalized critical nucleus size N c =N 0 c and the normalized interfacial energy γ/γ 0 on the degree of preordering (i.e., local bond order fluctuation) that is controlled by Δt MD . Here N 0 c and γ 0 are the corresponding values in the normal liquid state (without OKS). With shorter Δt MD , there is less preordering, leading to the increase of γ effectively. The dashed curve is an eye guide.
In order to further illustrate the crucial role of preordering in crystal growth, we use a two-step method to tune the structural ordering. As shown in Fig. 3b, we choose initial states with different crystal sizes from the simulations of Δt MD = 32 ps in Fig. 3a, and then we decrease Δt MD to 10 ps to kill the preordering more effectively. Obviously, more severe suppression of bcc-like preordering leads to the slower growth of the crystal. For the same Δt MD , the impact of order killing is weaker for larger crystals. In other words, in the late stage of crystal growth, we need to use shorter Δt MD to slow down the growth because the preordering is prevailing around the large crystal seed (see Fig. 1). This result indicates the curvature dependence of the interface preordering, i.e., the liquid-crystal interfacial energy. This demonstrates that the interface preordering is of paramount importance also during crystal growth. If there is no preordering wetting the liquid-crystal interface, the order-parameter profile across the interface between the crystal and the liquid becomes very sharp. Such a significant structural difference across the crystal growth front, i.e., the large interface tension, makes the crystal growth into the liquid more difficult (see the additional simulations on Si below). This is intuitively natural since a more significant structural change at the growth front is required for crystals to grow. CNT has overlooked this fact.

Critical nucleus size under the control of structural ordering
Inspired by the time dependence of N xtal on Δt MD during crystal growth, we measure the critical nucleus size N c under the conditions when preordering is suppressed by OKS (see Methods, section 'Order killing strategy'). Figure 3c illustrates the dependence of the N c normalized by N 0 c on the degree of preordering (i.e., local bond orientational order fluctuation) controlled by Δt MD . Intriguingly, with decreasing Δt MD , i.e., fewer preordering, N c =N 0 c increases more steeply, indicating the ascending critical nucleation barrier. Note that the value N c =N 0 c should plateau to 1.0 at large enough Δt MD eventually. The required time to reach the plateau will depend on the crystallisation time of the supercooled liquid. When Δt MD reaches the time when the crystal nucleus size exceeds the critical one, which we call the crystal nucleation time here, the order-killing will not affect the crystal nucleation process anymore. We estimated the average crystal nucleation time for NiAl at T nose as 7.38 ± 5.61 ns from 10 independent simulations (see Supplementary Fig. 4). This is much longer than the largest Δt MD value (0.2 ns) used in Fig. 3c, leading to the unsaturation of N/N 0 . Since the driving force for crystallisation Δμ is roughly unaltered (see Methods, section 'Crystallisation driving force'), the rise of N c should mainly originate from the increasing interfacial energy γ. In fact, CNT predicts the critical nucleus size as R = 2γ/(ρ s |Δμ|), where ρ s is the number density of the solid 1 . By assuming a spherical shape of the nucleus as in CNT, even though it is not necessarily to be the case, the change of γ normalised by the interface tension of the normal liquid state γ 0 scales as the cube root of N c =N 0 c . As a consequence, γ/γ 0 is estimated to increase from approximately 1.30 for Δt MD = 0.2 ns to around 1.82 when Δt MD decreases to 25 ps (Fig. 3c). Therefore, the way by which preordering tailors the crystallisation kinetics is by adjusting γ. Both crystal nucleation and growth proceed via the non-classical pathway mediated by the local bond order because the emergence of preordering at the liquid-crystal interface can effectively reduce γ by wetting the crystal and lower the crystallisation barrier.
The correlation between γ and Δt MD shows the temperature (T) dependent nature of the interfacial energy, different from CNT in which it is assumed as constant. This is because CNT assumes that a liquid has a completely disordered structure. γ(T) can either decrease or increase with supercooling until the glass transition intervenes 39 . This feature depends on the local structural ordering characteristics. In systems where crystal-like preordering are prevalent, lowering T will enhance these crystal-friendly structures 16 and reduce γ(T). This picture is consistent with the diffuse interface theory and has been confirmed by the studies on liquid metals [40][41][42] . The role of icosahedra in quasicrystal formation is similar to that of crystal-like preordering in ordinary crystallisation 28,29,37,43 . Such examples have been experimentally observed for ZrCu-based and Al-based bulk metallic glasses 44,45 . On the contrary, if local icosahedral order is dominant, but not strong enough for quasicrystal formation, and suppresses crystal-like preordering, γ(T) for crystal formation would increase and helps glass formation.

Phenomenological modifications to the classical theory
We have shown that local structural orderings are crucial in crystallisation, but which have been overlooked in the classical theories. In the following, we make some essential modifications of CNT to incorporate the effect of such orderings. The "homogeneous" nucleation frequency of crystals in supercooled liquids is rewritten phenomenologically as where C I is a constant prefactor, D is the translational diffusion coefficient normal to the crystal surface, k B is Boltzmann constant and α is a factor accounting for the shape of the crystal nucleus. Both γ and Δμ should have non-trivial temperature dependences. More importantly, γ should be closely related to the specific structural orders formed in a supercooled liquid and is determined by competing ordering effects 38,46 .
Regarding crystal growth, only two factors have been mainly considered traditionally: 1) the effective diffusion coefficient, which represents the atomic mobility towards the interface; 2) the driving force Δμ, which dictates the willingness of an atom to join the crystal. However, based on the above findings of the critical role of structural orderings in crystallisation, we phenomenologically introduce a new factor P(γ(T)) to account for the microscopic details of the interfacial structural ordering. This factor measures how easily an atom finds a proper lattice position on the crystal surface. The rate of crystal growth is accordingly modified as where ℓ and λ are respectively the inter-planar spacing and the displacement length based on the classical Wilson-Frenkel theory 1,26,27 . U WF (T) represents the Wilson-Frenkel theory without considering any effect from γ. We note that U WF (T) does not simply represent the crystal growth rate in the absence of preordering. It only considers the diffusion rate and thermodynamic driving force for crystallisation, neglecting the microscopic details of system-specific liquid-crystal interface properties. U WF (T) can strongly overestimate the crystal growth rate if the topological and chemical frustration effects are serious in a system. For example, in binary and multicomponent systems with large structural and chemical (compositional) gradients near the crystal surface, even without preordering, the crystal growth could be somewhat slower than the WF prediction. This may explain the discrepancy between experiments and theories previously reported 2,4 . To solve this issue, we add a new factor P(γ(T)) incorporating the effect of the liquid structural ordering of a specific system on the kinetic and/or thermodynamic factors. In some sense, P(γ(T)) measures the wettability of the interface preordering to the crystal, even though its absolute value does not directly give the degree of wettability. We stress that it contains the kinetic effects. A large P(γ(T)) above unity indicates that crystal-like preorders wet the crystal. Thus, the WF theory underestimates the crystal growth rate. This scenario might explain the ultrafast crystal growth in pure metals far exceeding the prediction of the classical theories 21 . However, the meaning of P(γ(T)) < 1.0 should depend on the actual situation. This may be related to the temperature dependence of γ in a system 39 . If γ(T) decreases with increasing degree of supercooling, more crystal-like preordering could form when the temperature is lowered. Thus, the crystal phase should be more wetted. So, even when P(γ(T)) < 1.0, there may be wetting of crystals by preorders. However, if there is no preordering during cooling, P(γ(T)) < 1.0 tells no wettability. This is because if the wetting effect exists, the crystal growth should be accelerated.
To verify our idea, we directly calculate the variables in equation (2) to estimate the WF prediction and the actual crystal growth rate from MD simulations (see Methods, section 'Crystal growth rate measurement'). The properties of the eight systems studied are listed in Table 1. The systems cover from single-component systems to binaries, metals to metalloids, isotropic interactions to anisotropic interactions, and three-body interactions to many-body interactions. They also have different equilibrium crystal structures and different types of locally favoured structures. Thus, the systems have enough variety for a general discussion valid to a broad range of materials. We show in Fig. 4 the temperature dependence of the measured crystal growth rate of each system (solid symbols), together with the values of the WF prediction (red dashed line). The WF theory fails seriously in describing the measured crystal growth rate. It underestimates U of monoatomic metals (Ta, Zr, Cu and Pt) and overestimates the ones of binaries and tetrahedral systems (CuZr, NiAl, Si, and water). The former case has also been verified in several fcc-forming metals 21 . In addition, besides water and Si, the temperature dependence of the WF prediction is distinctly different from the measured values. We show the temperature dependence of the correction factor P(γ(T)), which measures the discrepancy between the WF prediction and the measured U value, in Fig. 5. This result suggests that γ depends on the supercooling degree, as inferred from Fig. 3c. Increasing P(γ(T)) with lowering temperature implies a decreasing γ with supercooling. In the metallic systems, crystal-like bond orientational ordering formed in the supercooled state, as revealed in NiAl above, should help crystal growth. This effect is especially prominent in monoatomic liquid metals, which do not suffer from chemical (compositional) frustrations 38,[47][48][49][50] . For Si and water, the local tetrahedral order will facilitate the growth of the diamond cubic crystal phase. These wetting effects originate from the compatible local symmetry between the preordering and the crystal phase.
The temperature dependence of P(γ(T)) in Fig. 5 also shows some interesting features. On the one hand, there is a broad range of P(γ(T)) for these systems, demonstrating different characteristics of liquid structural orderings when assisting crystal growth. On the other hand, the dependence of P(γ(T)) on the degree of supercooling is significantly different among the systems. P(γ(T)) increases when more supercooled for monoatomic metals. The binary NiAl behaves somewhat similar to metals but exhibits a weaker increase of P(γ(T)) with decreasing T than metals. However, for CuZr, Si, and water, P(γ(T)) keeps almost constant at various degrees of undercooling. We analyse the structural orderings in the equilibrated supercooled states and at the crystal-growth front for Zr, Ta, NiAl, CuZr, and Si to rationalise these observations. First of all, the degree of the spatial correlation of crystal-like preordering is important in crystal growth. A stronger correlation would help crystal growth more strongly. We compare the spatial correlation of the bond orientational order parameter (Q 6 and Q 12 ) in Supplementary Fig. 11 and find a substantial increase in the spatial correlation of the preordering upon cooling for metals. Such correlation is weaker for NiAl and even weaker and almost temperature independent for CuZr. In Si, the spatial correlation of the order parameter is nearly absent. Secondly, we show the spatial distribution of various structural orderings at 0.7T m , below which Zr is hard to equilibrate without crystallisation, in Fig. 6. As we can see, in the supercooled state, there are lots of ico-like orderings in Ta, while bcc-like preordering is dominant in Zr (Fig. 6a, b). Both show rich preordering at the liquidcrystal interface. However, there are much more ico-like orders adjacent to crystal-like preorders in Ta than Zr. The ico-like ordering should hinder the growth of crystal-like preorders into the liquid region. This feature explains the slower crystal growth rate (Fig. 4) and better glass-forming ability of Ta than Zr 51 . As shown in Fig. 6d, similar to Fig. 1b, some ico-like orders are coexisting with bcc-like preorders in the supercooled state. This situation is similar to CuZr, but which is more abundant of ico-like orders while the bcc-like preorders become rarer and more dispersed in space. This relationship between NiAl and CuZr is similarly found for the crystal-growth front structure (Fig. 6d, e,  middle panel). Even though the bcc-like preorders wet the interfacial crystal in CuZr, the large amount of ico-like orders adjacent to the  − 6/14), where c α is defined as the fraction of α-type atoms in the nearest neighbours of a central α-type particle (the white atoms are crystals). In the cases of crystal growth, only a part of the systems is shown for clarity. For the metallic systems, the pink, cyan, green, and blue colours represent bcc-like (preorder), ico-like, fcc-like, and hcp-like structural orderings, respectively. 14 or 12 nearest neighbours are shown for bcc-like or ico-or fcc-or hcp-like orders, respectively. In Si, the pink atoms have 10-11 connected neighbours, the cyan atoms have 6-9 connected neighbours, and the green atoms have 5 connected neighbours. Four neighbouring atoms are included for visualisation (i.e. forming a tetrahedral, for example). In all systems, the yellow atoms denote crystals. The sizes of particles have been adjusted for better visualisation. The normal liquid atoms are not shown. interface would impede crystal growth. Meanwhile, both NiAl and CuZr suffer from compositional frustration 38,47-50 (see also Supplementary  Fig. 12), which is even more severe in the latter (Fig. 6d, e, right panel). The bond strength between atoms should be very critical in determining the local composition adjustment rate 38 , which is directly linked to the crystal growth rate. Furthermore, the slower-growing speed, i.e., smaller P(γ(T)), of NiAl than the monoatomic metals (see Fig. 5) should originate from more ico-like orders and more substantial compositional frustration at the interface fronts in NiAl than these metals. This demonstrates the significant role of chemicalcomposition fluctuation in crystal growth, an intrinsic, unique feature of multicomponent (metallic) systems. This provides a physical explanation for the increase in the glass-forming ability of metallic alloys with the number of components. As for the tetrahedral system like Si (Fig. 6c), the preorders are much fewer and more independent from each other in space than monoatomic metals. Therefore, the liquid-crystal interface is much sharper. There is a steep gradient in the structural order parameter at the interface, without composition fluctuation, which makes P(γ(T)) small and independent of temperature change.
Thus, the difference of P(γ(T)) among these systems indicates that the crystal growth rate is determined by the degree of crystal-like preordering on the crystal surface and, additionally, the compositional fluctuation for the binary systems. All these microscopic features contribute to determining the liquid-crystal interfacial energy γ. How quantitatively each of these factors determines γ is an interesting topic for future study. These results strongly indicate the crucial role of interfacial preordering and interfacial energy in the crystal growth kinetics.

Discussion
Our findings suggest that the interface tension γ(T) could be the most crucial factor in governing the crystallisation kinetics and glassforming ability of various substances, including metallic glasses 20,38,46 . It means that the glass-forming ability is microscopically governed by the competing ordering effects 16 . Considering the case in NiAl, with Δt MD = 0.2 ns, the increase of γ(T nose ) only by a factor of 1.30 leads to a delay of crystallisation by nearly three orders of magnitude. As is well known, the crystallisation time at T nose is the shortest and directly determines the critical cooling rate. Therefore, suppressing the fluctuation of local bond orientational order friendly to the crystal improves the glass-forming ability significantly. Further enhancement of γ(T nose ) can even turn a poor glass-former into an excellent one. Our findings also indicate the different temperature dependences of γ in different materials, which may provide another crucial piece for improving CNT and better understanding the glass-forming ability in metallic glasses.
In conclusion, we have studied the crystallisation mechanism in a supercooled metallic liquid and found that crystal forms and grows via non-classical pathways. Both crystal nucleation and growth are mediated by local crystal-like bond orientational order fluctuation, seriously questioning the assumptions of the classical theories. By developing the new order-killing technique, we successfully manipulate the crystallisation kinetics to a large extent, even at the same thermal condition. The degree of liquid preordering controls crystallisation by tailoring the liquid-crystal interface's properties. This indicates the crucial role of interfacial energy in crystallisation and glass formation. We have shown that an essential modification to the classical theory is necessary, evaluated the correction factor, and discussed its meaning for various material types. Crystal-like preordering may influence the kinetic and/or thermodynamic factors. It is highly desirable to theoretically reveal the underlying physical mechanism of how preordering affects crystallisation kinetics. Our findings should deepen the fundamental understanding of the crystallisation mechanism in supercooled liquids and pave a new way to control crystallisation and glass formation.

MD simulations of NiAl crystallisation
We perform extensive MD simulations for NiAl by employing the many-body embedded atom method (EAM) potential 52 . All the simulations are run by using the open-source LAMMPS software 53 , in which periodic boundary conditions are kept in three directions and the time step for integral is 0.002 ps. Isobaric-isothermal ensemble (NPT) with zero external pressure is employed. We first create a bcc lattice with N = 8192 atoms and equilibrate it at a high temperature (2000 K). Then, we quench the system instantly to the desired temperature T relax (mainly the nose temperature T nose from the Temperature-Time-Transformation curve) for isothermal annealing and wait for crystallisation. The system size is designed so that there is only one critical nucleus during crystallisation. For NiAl, T nose ≊ 921K(~0.6T m ) (see Supplementary Fig. 1 and Ref. [54]). 8 independent simulations are generally carried out for all simulations, which give similar results. The simulation methods for other systems (see below) are essentially similar.

Order killing strategy
To effectively control the structural (or, more precisely, local bond orientational order) fluctuations, we develop OKS by periodically performing biased MD simulations and normal MD simulations. The biased MD simulations are performed by using PLUMNED 2 patched to LAMMPS 55 . During the isothermal annealing, we check the local bond orders at a period of Δt MD (normal MD simulation timescale in each round). Suppose preorders are detected by bond orientational order parameter Q 6 (see below). In that case, we perform MD simulations with a harmonic bias potential acting on these preorders to revert them to the disordered liquid state. This is realised by displacing the centre atoms. We emphasise that we do not bias crystallised atoms. If there are no preorders, we will run normal MD simulations. So Δt MD is key to control the degree of structural fluctuations. The shorter Δt MD is, the stronger the suppression of structural fluctuations will be.
To construct the bias potential, we first need to define a collective variable q 6v as the average q 6 (i) of the selected atoms. In the context of enhanced sampling of PLUMED 2, q 6 (i) is defined as where Y 6m is the sixth order spherical harmonics and is the switching function to determine whether atom j is the nearest neighbor of atom i. We choose r 0 to be 3.5 Å for better efficiency. The bias potential takes the form of in which we set k = 10, s = 0.1 and a = 0. The restraining potential (UPPER_WALLS) starts acting on the preorderings when their q 6v is above 0 and tries to decrease it back to 0. The effective time of the bias potential is set as 10 4 MD steps, which is large enough for energy dissipation in our case. In this way, we are able to suppress the structural fluctuations in a controllable manner. OKS is a little bit computational costly because we need to characterize the local structure frequently and calculate the collective variable every MD step during the biased simulations. When we study crystal growth or measure the critical nucleus size under a bias, a crystal seed of a specific size is first generated in the supercooled liquid by normal MD simulations in spontaneous crystallisation. We then use OKS during the seeding simulations. The critical nucleus is identified when the probability of the crystal seed growing or dissolving becomes equal, i.e., when the crystal seed becomes stable. As described below, we identify the crystallised atoms through the bond orientational order methods for the crystal nucleation and growth processes.
The effective force to kill preorder from the implementation of V bias may increase the temperature slightly in the liquid state. However, since the bias force is relatively weak, we confirm that such temperature effects on the crystal nucleation and growth are very weak (see Supplementary Figs. 5 and 6). Similar results are also observed for the fluctuations of pressure and number density (see Supplementary  Fig. 7). Unlike the conventional biasing method, where one can clearly separate the free energy part and the kinetic factor part, our strategy involves the kinetic factor in the biasing. Nevertheless, we note that the influence of the bias potential on the atomic motion is negligible. As shown in Supplementary Fig. 8, the atom position is only perturbed inside its cage by the biasing. Thus, the periodic perturbation does not affect the rate factor, i.e., particle diffusion. Therefore, we can conclude that OKS only affects the thermodynamic factor through perturbating crystal-like orders in the liquid without affecting the kinetic factor.

Structural characterisation
We use bond orientaional order parameters to define the local structures in supercooled metallic liquids 56 . The neighbours of each atom is defined by the radical Voronoi tessellation 57 . Firstly a complex vector for atom i is computed as where N i is the number of the nearest neighbours of atom i, − l ≤ m ≤ l, and Y lm is the spherical harmonics. f j is the fraction of Voronoi polyhedral face area between i and j over the overall face area. The coarsegrained order parameter Q l as in which Q lm (i) is the average of q lm (i) over the first coordination shell, including atom i itself. To detect the crystallised atoms in the metallic systems (Ta, Zr, Pt, Cu, NiAl and CuZr, see below), we define an order parameter In this calculation, we weighted the spherical harmonics of each atom from its nearest neighbours by their Voronoi face area, respectively 58 . The bond between i and j is treated as crystalline if S 6 (i, j) > 0.7, and if the number of crystalline bonds exceeds 10 we treat i as crystallised. As for NiAl, specifically, since the primary crystalline phase is B2 (bcc-like), we find there are no fcc-like or hcp-like structural orderings during crystallisation, which is different from the colloidal systems 7 . Therefore, the atoms with Q 6 larger than 0.25 but with the number of crystalline bonds less than 10 are considered preorders (with bcc-like symmetry). We use the parameter w 6 defined as to identify the icosahedral center (w 6 < − 0.023) 59 . The term in parentheses is the Wigner 3-j symbol. The different crystal-like structural orderings (bcc-like, fcc-like, hcp-like) are characterised by the order parameters in which l ∈ {4, 6}. A particle belongs to bcc-like crystal is indicated by W 6 > 0, while for fcc-like crystalline atoms (W 6 < 0, W 4 < 0) and for hcplike crystalline atoms (W 6 < 0, W 4 > 0). Similarly, for tetrahedral systems (Si and water), we define S 12 (i, j) as the following to characterise the crystallised atoms 46 : We take the nearest 16 atoms as the particle neighbours in these calculations. An atom is treated as crystallised if it has more than 12 crystalline bonds (connected neighbours) as S 12 (i, j) > 0.75 for a pair.
To characterise the properties of the structural ordering in space, we calculate the spatial correlation function based on Q l (l = 6, 12) as Upon selecting preorders from the liquid state to kill, we choose the bcc-like preorders by Q 6 > 0.25 in which the spherical harmonics are not weighted by the Voronoi face area. This simple method is effective because the fcc-like and hcp-like preorders are negligible (see Fig. 1). It is beneficial in saving computational time. On the other hand, the local icosahedral orderings can be automatically excluded (see Supplementary Fig. 9).

MD simulations on additional systems
To study the crystal growth kinetics, we perform additional simulations on eight different systems: Ta, Zr, Cu, Pt, CuZr, NiAl, Si and water. The first six systems are modeled by the empirical many-body EAM potentials 51,52,60,61 , while the interatomic interactions in the latter two systems are described by the three-body Stillinger-Weber potentials 62,63 . All the simulations are run by using the open-source LAMMPS software 53 , in which periodic boundary conditions are kept in three directions, and the time step for integral is 5 fs for water and 0.002 ps for the others. The specific simulation strategies are designed for different measurements (see below). Generally, five independent simulations are performed for ensemble average.

Crystallisation driving force
In our simulations, the thermodynamic driving force for crystallisation Δμ(T) is estimated empirically via where ΔH m is the enthalpy of fusion per atom at the melting temperature T m . To measure ΔH m (T m ) we create equilibrium solid (crystal) and liquid at T m respectively and calculate the enthalpy difference directly. We first create a cubic cell based on the specific material crystal structure (see Table 1). We generate 18 × 18 × 18 unit cells for bcc and B2 types, 14 × 14 × 14 unit cells for fcc type, and 11 × 11 × 11 unit cells for dc type. So roughly N~11000 atoms are considered for each system. NPT (constant number N, constant pressure P, constant temperature T) emsemble is used for all the simulations. The simulation procedure consists of the following steps: (1) relax the initial configuration at a temperature much lower than T m ; (2) heat the crystal to T m ; (3) equilibrate the crystal at T m ; (4) heat the crystal to a temperature much above T m ; (5) equilibrate the melt at that high temperature; (6) cool the liquid to T m ; (7) equilibrate the liquid at T m . The simulation time of step (5) is 4.0 ns, and it is 1.0 ns for the other steps. The enthalpy of the crystal and liquid is calculated from steps (3) and (7), respectively.

Crystal growth rate measurement
We construct a planar liquid-crystal interface (see Supplementary  Fig. 10 as an example) to calculate the crystal growth rate U(T) at various temperatures for the eight systems. A rectangular cell is designed for each system based on the material crystal structure (see Table 1). Specifically, we create 36 × 12 × 12 unit cells for bcc and B2 types, 30 × 10 × 10 unit cells for fcc type, and 24 × 8 × 8 unit cells for dc type. The liquid-crystal coexistent configuration is then generated by pinning part of the crystal while melting the others after equilibrating the initial crystal configuration. After relaxing the coexistent configuration at T m , the system will be annealed at the target temperature to monitor crystal growth. In the annealing process, NPT ensemble is used, but only the pressure along the elongated direction is controlled at 0, so only the box length along that direction can be adjusted for crystal growth and melting. The crystallisation process is analysed by using S 6 (i, j) or S 12 (i, j) as described above. The crystal growth rate is estimated by the crystal growth distance normal to the initial interface and the corresponding time elapsed.

Diffusion coefficient
The diffusion coefficient D(T) is calculated from the mean-squared displacements To measure D(T), we equilibrate the liquids at each temperature and then output atomic trajectories to calculate the mean-squared displacements (see examples in Supplementary Fig. 13). In more detail, we first equilibrate the liquid at a temperature much higher than T m and then bring the configuration to the target temperature for further equilibration. The ensemble will then be switched from NPT to NVT (constant number N, constant volume V, constant temperature T) for production runs. We use N = 11664 atoms for these calculations. The configurations from these simulations are also used to analyse the structural orderings in the equilibrated supercooled liquids.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/.